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Abstract 

Modern dynamical systems theory has previously had little to say 
about finite difference and finite element approximations of partial 
differential equations (pdes) Q. However, recently I have shown one 
way that centre manifold theory may be used to create and support 
the spatial discretisation of pdes such as Burgers' equation Q and 
the Kuramoto-Sivashinsky equation |J. In this paper the geometric 
view of a centre manifold is used to provide correct initial conditions 
for numerical discretisations [Qj . The derived projection of initial con- 
ditions follows from the physical processes expressed in the pdes and 
so is appropriately conservative. This rational approach increases the 
accuracy of forecasts made with finite difference models. 
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1 Introduction 

Consider the equations for some physical field u(x,t) evolving in space-time 
that we wish to model numerically. Imagine a given initial field uq(x) and 
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a finite difference model written in terms of %(t) = u(xj,t) for equi-spaced 
grid points Xj — jh say; for example, in §||] for Burgers' equation (§) we find 

in terms of the central difference operator Suj = Wj+1/2 — u j-i/2 and central 
mean operator fiUj = (uj+i/z+Uj-ip) /2 . One might expect that the correct 
initial condition for this discretisation is simply to project the initial field 
uq(x) onto the finite dimensional space of the model by setting the initial 
discretisation values to the value of the initial field at the grid: Uj(0) = 
uo(xj). But if the initial field is localised away from any grid point then 
physically we know to distribute the initial field among nearby grid points. 
I use dynamical arguments to show that the correct initial condition is, to 
leading order, the correctly conservative element average 

I rXj+h/2 

Uj(0) « — / uq(x) dx . (2) 

h Jxj-h/2 

This formula, and higher order corrections that involve neighbouring ele- 
ments, are derived systematically herein. For a numerical model, this is the 
first time a dynamical rationale has been used to provide initial conditions. 

Such projection of initial fields onto the discretisation is supported by 
centre manifold theory H e.g.]: the Relevance Theorem asserts that each 
of the nearby solutions of the governing pde exponentially quickly in time 
approach a solution of the numerical model; this holds even for finite grid 
spacing h. The algebraic techniques developed by Roberts M, based upon 
analysing with the aid of computer algebra the adjoint of a linearisation of the 
pde, determines the initial condition for the discretisation so that we ensure 
the finite difference model faithfully tracks the correct particular solution of 
the pde. 



2 Burgers' equation is discretised with centre 
manifold theory 

Consider the dynamics of Burgers' equation 

u t + auu x = u xx (3) 

as a prototype advection-diffusion equation. Roberts @ first constructed 
finite difference approximations to the spatial derivatives using centre man- 
ifold theory to ensure nonlinear, subgrid-scale processes were systematically 
modelled. We summarise the approach in this section. 
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Divide the spatial domain / into a number, say m, of elements of equi- 
size h. We analyse the dynamics of the elements away from any physical 
boundary to derive a discretisation for the interior of the domain. Artifi- 
cially crafted internal boundary conditions (ibc's) between the elements are 
introduced: 



u x — \au 2 







(1 - j)h(u x - \av?) = 7 [u] 



(4) 



where [ ] denote the jump across each internal boundary, — denotes the 
average value from the two sides of the boundary, and distinct from earlier 
work P| these ibc's are expressed in terms of the flux q = —u x + |aw 2 . See 
that when 7 = the right-hand side of the second ibc disappears so that 
the two conditions then completely insulate an element from its neighbours. 
Whereas when 7 = 1, the left-hand side disappears and the two conditions 
ensure sufficient continuity of the physical field to recover Burgers' dynamics 
throughout the domain. 

The centre manifold and the evolution thereon is straightforwardly con- 
structed using the computer algebra algorithm described in [TJ, |6|. Here we 
find the subgrid field in the jth element is 

u(x, t) = uj + ^ah^u 2 + ^a 2 h 2 ^ 2 u 3 + 7 £fi5uj + \i 2 5 2 Uj 
+ a/ry —^(uj5 2 Uj + 5 2 u 2 + 4u 2 ) + |£ 2 (2uj/i5 Uj — fi5u 2 ) + ^ 2 Uj8 2 u 3 



+ a 2 h 2 ^ jQ^(ujfi5 u 2 + fx5 



/<„- 



C(3u 2 5 2 u 



2ujS 2 Uj 



S 2 ^ + 8u 3 A 



+ \e{2u 2 ii8ui - u^5u 2 ) + U 4 u 2 S\] + O (7 



2 a 3 



(5) 



where £ = (x — Xj)/h ranges over [—1/2, 1/2] . The evolution on this centre 
manifold, when evaluated at 7 = 1 to restore continuity, forms the finite 
difference model (|l|) for Burgers' equation: see that the first three terms in ([I]) 
form a standard discretisation of each term but now appearing automatically 
from the discretisation when mediated by the flux form (^) of the ibc's; 
whereas the last term gives O (a 2 ) corrections to account for interactions 
between the nonlinear advection and the diffusive dissipation. Such nonlinear 
modifications of standard discretisations can be extremely effective 0. 

To find the correct initial condition, Uj(0), for numerical models such 
as (HD corresponding to any given field uq(x), we follow the procedure de- 
scribed in [|J. The aim is to determine projection vectors z 
those shown in Figure [I], so that 

1 

hJi 

as the inner product. Now the dynamics linearised about the nonlinear centre 



n x), such as 



(zj, Uo(x) — v(u(0), x)) — using (z,u) — — [ zudx 

h J 1 



(6) 



3 



T 



1 " 




-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 

^(x-x.)/h 

Figure 1: leading two orders of approximation to the projection vectors Zj(x) 
for purely diffusive dynamics, errors: O (7), dashed; O (7 2 ), solid. 



manifold, u = v(u,x), is governed by the operator 

J = d 2 x - av x - avd x , 

with ibc's linearised about (Q) of 

[u x } = 0, (l-j)h(u x -avu) = j[u], (7) 

Then in the above inner product the adjoint of J is 

J*z = d x z + avd x z , such that [z x ] = , (1 — 7)/^ = 7 [z] . (8) 

To find the projection vectors Zj(x) we start with the leading approximation 
Zj(x) ~ Xj{ x ) corresponding to @ and plotted in Figure |I], where Xj( x ) 
denotes the characteristic function that is 1 in the jth element and otherwise 
is 0. Then successive corrections are sought by iteration to ultimately satisfy 
the appropriate version of the equations derived in ||: defining the dual 
operator Vz = ff + J^z we must solve 

Vzj-Y t (Vz ji e i )zi = Q t (9) 

i 

subject to the ibc's in and the normalisation condition 

(zj,ei) = 6ij , (10) 

where ej = dv/duj is the tangent vector of the centre manifold. We seek 
solutions in a power series in 7 to errors O (7^ corresponding to the finite 
difference approximation of stencil width 2£—l. A computer algebra program 
available from the author does all the necessary algebra. 
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3 Project onto Burgers' discretisation 

In this section we solve to quantities with errors (9(a 3 ,7 2 ): the finite dif- 
ference model for Burgers' equation is then (|l[); and the corresponding cen- 
tre manifold over the whole domain is given by (^). Calculating to errors 
O (a 3 ,7 2 ) the projection onto the numerical model must be orthogonal to 



h 2 a 2 9 

1 " ~w~ > 



Xj 



+ 7 [(I - e) X 3 + + + K 2 ) Xi-i + K 2 ) Xj+i 



+ 



(-(12f - 16£ 3 )u, + Uj+1 - Wy_i) Xi 



+ 



ft, 07 
"48" 

+ - (3 - 6e + 8e 3 )^-i j 

+ + (3 + 6£ - 8Z 3 )u J+1 ) Xj+i 

h 2 a 2 ^ 
384 



(8(1 + 3£> 2 + A Uj (u j+1 + Uj-x) + 2{u 2 +l + u 2 ^)) Xj 



+ (3u 2 + 4ujUj-_i - (5 + 12£ 2 + le^ 3 )^) Xi-i 



+ 3« 2 + 4 W+1 - (5 + 12£ 2 - 16£> 2 



(11) 



Higher order expressions may be straightforwardly computed by computer 
algebra. I conclude by further interpreting the physical effects incorporated 
in the projection defined by the above Zj. 



3.1 Linear diffusion 

Set a = in this subsection to analyse the linear diffusion equation u t = u xx . 
Then the projection vector (|TT|) , evaluated at 7 = 1 to recover the physically 
relevant plotted in Figure [I], is 

% « (I - £ 2 ) X, + + \i + l^ 2 ) tt-i + ~U+ ^) X j+ i ■ (12) 

To find the correct initial condition using this in note that in these linear 
diffusion dynamics (zj, v(u, x)) = Uj by the normalisation (|I~OD; thus here 
Uj(0) = (zj,uq(x)) . For example, see that a point release in the fcth element, 
uq(x) = S(x — Xk — hrj), requires the slightly distributed initial condition 

ft^-(O) = (I - V 2 ) 8 ktj + (-i -\ V + \rf) + (-£ + \r) + W) • 

(13) 

Such a specific initial condition corresponds via ([5]) to a field on the centre 
manifold as shown in Figure ^for the three cases 77 = 0, 1/4 and 1/2. 
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Figure 2: initial fields u = v(u(0),x) corresponding to a unit-mass point 
release at: £ = 0, dot-dash; £ = 1/4, dashed; £ = 1/2, solid. 



See that these initial conditions ensure that the first moment of the nu- 
merical solution is correct for all time: in the numerical model ([3]) the first 
moment is constant in time so it is enough to check that the first moment is 
correct in the initial conditions. Define (u) = (1, u) , then for all time (u) = 1 
both in the model and in the exact solutions. The first moment in the exact 
solution is its initial value mi = ((x — Xk)u (x)) = hrj; from ( JTJ ) the first 
moment in the numerical model is the same 

m i — ({x — x k) v i u {fy, x)) = hrj . 

However, the second moment rri2 = ({x — Xk) 2 u) has O (h 2 ) errors: it evolves 
in time at the correct rate drri2/dt = 2, but the initial value is h 2 (r] 2 — 1/6) 
instead of 0. Determining the projection of initial conditions to higher orders 
in the coupling parameter 7 obtains such higher order moments correctly. 
Note that the rational approach adopted here does better than the usually 
chosen initial conditions which incur O (h) errors. 

3.2 Nonlinear dynamics 

Consider the O (a) terms from (|ll|) that modify (|T2"D , namely 

ha 



2 = + 



3 



48 



((-12f + 16£ 3 )n, + u j+1 - u^) Xi 
+ (+Uj + (-3 + 6£ - S^ 3 )^) Xj-i 
+ (- Uj + (+3 + 6£ - 8£ 3 )u J+1 ) Xj+i ■ 
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Realise that the leading order effect of including these terms is to modify the 
initial condition by (z'^u^x)^ . For example, if the initial field is approxi- 
mately constant, u (x) « U, then 

z 'i = ^ ~ 8 ^)(^'+i - 2 *i + Xj-i) + 2(Xi+i " Xi-0] ; 

that the coefficients of the characteristic functions Xk sum to zero reflects 
that the the projection conserves the field u. More specifically, if uq(x) 
is U except for a symmetric bump in the kth element, then as well as the 
direct symmetric distribution identified for linear diffusion, the component in 
Xj+i — Xj-i causes Uk-i{0) to increase and iz fc+1 (0) to decrease by an amount 
proportional to Ua reflecting that the self advection of the bump is not as 
great as that induced by assigning the mass of the bump solely to Uk(0). 
Conversely, for an antisymmetric perturbation in the kth element, positive 
to the left of Xk, the component in Xj+i ~ %Xj + Xj-i increases Uk(0) and 
decreases M fc±1 (0) in proportion to Ua to reflect the increased delay in u 
advecting out of the kth element because more of it is further to the left 
initially. The O (a 2 ) terms in (|TT| ) reflect more subtle physical processes. 



4 Conclusion 

Based upon the method of analysis and the discussion in the previous sec- 
tions, we deduce that this centre manifold approach to finding correct initial 
conditions for finite difference models accounts for subgrid scale processes 
that occur as initial transients decay. No other method does this. 

Extensions of this approach to higher spatial dimensions is straightfor- 
ward. For example, consider the class of diffusive pde's 

^ = V 2 u + f(u,Vu), 

where / represents nonlinear reaction or advection effects. After tessellating 
space into finite elements — using ibc's of the form (cf (|])) 

[q n ] = and (1 - j)hcfc = j[u] 

where q n is the flux of u normal to the internal boundary and h is a size 
of the element — the fundamental problem in constructing a model is simply 
to solve Poisson's equation with forced Neumann boundary conditions on 
each element. The adjoint of this problem lies at the heart of the dual 
for determining initial conditions of the approximation. Although these sub- 
grid problem may itself need to be done numerically, in the simplest case of a 
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regular tessellation it need only be done once for each term in the model, just 
like the computation of the interaction terms in a traditional finite element 
approximation. 

In the case where there are variations in the size or shape of the ele- 
ments of the discretisation, one would build formulae for the approximation 
parametrised by the shapes of the element and those neighbouring elements 
to which it is coupled by the ibc's. The algebraic detail becomes more com- 
plicated but the principles are the same. 
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